<!-- saved from url=(0014)about:internet -->
<html>
<body>
<script src="resources/eml_report_loadable_data.js"></script>
<div class="dead">
<pre class="code">
<span class="srcline"><span class="lineno"><a href="1,1" id="srcline1">  1</a></span><span class="line">function fused = Ying_2017_CAIP(I, mu, k, a, b) <span class="comment">% camera a, b</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,2" id="srcline2">  2</a></span><span class="line"><span class="comment">%%</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,3" id="srcline3">  3</a></span><span class="line"><span class="comment">% @inproceedings{ying2017new,</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,4" id="srcline4">  4</a></span><span class="line"><span class="comment">%   title={A New Image Contrast Enhancement Algorithm Using Exposure Fusion Framework},</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,5" id="srcline5">  5</a></span><span class="line"><span class="comment">%   author={Ying, Zhenqiang and Li, Ge and Ren, Yurui and Wang, Ronggang and Wang, Wenmin},</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,6" id="srcline6">  6</a></span><span class="line"><span class="comment">%   booktitle={International Conference on Computer Analysis of Images and Patterns},</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,7" id="srcline7">  7</a></span><span class="line"><span class="comment">%   pages={36--46},</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,8" id="srcline8">  8</a></span><span class="line"><span class="comment">%   year={2017},</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,9" id="srcline9">  9</a></span><span class="line"><span class="comment">%   organization={Springer}</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,10" id="srcline10"> 10</a></span><span class="line"><span class="comment">% }</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,11" id="srcline11"> 11</a></span><span class="line"><span class="comment">%</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,12" id="srcline12"> 12</a></span><span class="line"><span class="comment">% Please feel free to contact me (yingzhenqiang-at-gmail-dot-com) if you</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,13" id="srcline13"> 13</a></span><span class="line"><span class="comment">% have any questions or concerns.</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,14" id="srcline14"> 14</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,15" id="srcline15"> 15</a></span><span class="line">if  ~exist( 'mu', 'var' )</span></span>
<span class="srcline"><span class="lineno"><a href="1,16" id="srcline16"> 16</a></span><span class="line">    mu = 0.5;</span></span>
<span class="srcline"><span class="lineno"><a href="1,17" id="srcline17"> 17</a></span><span class="line">end</span></span>
<span class="srcline"><span class="lineno"><a href="1,18" id="srcline18"> 18</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,19" id="srcline19"> 19</a></span><span class="line">if ~exist( 'a', 'var' )</span></span>
<span class="srcline"><span class="lineno"><a href="1,20" id="srcline20"> 20</a></span><span class="line">    a = -0.3293;</span></span>
<span class="srcline"><span class="lineno"><a href="1,21" id="srcline21"> 21</a></span><span class="line">end</span></span>
<span class="srcline"><span class="lineno"><a href="1,22" id="srcline22"> 22</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,23" id="srcline23"> 23</a></span><span class="line">if ~exist( 'b', 'var' )</span></span>
<span class="srcline"><span class="lineno"><a href="1,24" id="srcline24"> 24</a></span><span class="line">    b = 1.1258;</span></span>
<span class="srcline"><span class="lineno"><a href="1,25" id="srcline25"> 25</a></span><span class="line">end</span></span>
<span class="srcline"><span class="lineno"><a href="1,26" id="srcline26"> 26</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,27" id="srcline27"> 27</a></span><span class="line">if ~isfloat(I)</span></span>
<span class="srcline"><span class="lineno"><a href="1,28" id="srcline28"> 28</a></span><span class="line">    I = im2double( I );</span></span>
<span class="srcline"><span class="lineno"><a href="1,29" id="srcline29"> 29</a></span><span class="line">end</span></span>
<span class="srcline"><span class="lineno"><a href="1,30" id="srcline30"> 30</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,31" id="srcline31"> 31</a></span><span class="line">lambda = 0.5;</span></span>
<span class="srcline"><span class="lineno"><a href="1,32" id="srcline32"> 32</a></span><span class="line">sigma = 5;</span></span>
<span class="srcline"><span class="lineno"><a href="1,33" id="srcline33"> 33</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,34" id="srcline34"> 34</a></span><span class="line"><span class="comment">%% t: scene illumination map</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,35" id="srcline35"> 35</a></span><span class="line">t_b = max( I, [], 3 ); <span class="comment">% also work for single-channel image</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,36" id="srcline36"> 36</a></span><span class="line">t_our =  imresize( tsmooth( imresize( t_b, 0.5 ), lambda, sigma ), size( t_b ) );</span></span>
<span class="srcline"><span class="lineno"><a href="1,37" id="srcline37"> 37</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,38" id="srcline38"> 38</a></span><span class="line"><span class="comment">%% k: exposure ratio</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,39" id="srcline39"> 39</a></span><span class="line">if  ~exist( 'k', 'var' ) || isempty(k)</span></span>
<span class="srcline"><span class="lineno"><a href="1,40" id="srcline40"> 40</a></span><span class="line">    isBad = t_our &lt; 0.5;</span></span>
<span class="srcline"><span class="lineno"><a href="1,41" id="srcline41"> 41</a></span><span class="line">    J = maxEntropyEnhance(I, isBad);</span></span>
<span class="srcline"><span class="lineno"><a href="1,42" id="srcline42"> 42</a></span><span class="line">else</span></span>
<span class="srcline"><span class="lineno"><a href="1,43" id="srcline43"> 43</a></span><span class="line">    J = applyK(I, k, a, b); <span class="comment">%k</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,44" id="srcline44"> 44</a></span><span class="line">    J = min(J, 1); <span class="comment">% fix overflow</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,45" id="srcline45"> 45</a></span><span class="line">end</span></span>
<span class="srcline"><span class="lineno"><a href="1,46" id="srcline46"> 46</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,47" id="srcline47"> 47</a></span><span class="line"><span class="comment">%% W: Weight Matrix </span></span></span>
<span class="srcline"><span class="lineno"><a href="1,48" id="srcline48"> 48</a></span><span class="line">t = repmat(t_our, [1 1 size(I,3)]);</span></span>
<span class="srcline"><span class="lineno"><a href="1,49" id="srcline49"> 49</a></span><span class="line">W = t.^mu;</span></span>
<span class="srcline"><span class="lineno"><a href="1,50" id="srcline50"> 50</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,51" id="srcline51"> 51</a></span><span class="line">I2 = I.*W;</span></span>
<span class="srcline"><span class="lineno"><a href="1,52" id="srcline52"> 52</a></span><span class="line">J2 = J.*(1-W);</span></span>
<span class="srcline"><span class="lineno"><a href="1,53" id="srcline53"> 53</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,54" id="srcline54"> 54</a></span><span class="line">fused = I2 + J2;</span></span>
<span class="srcline"><span class="lineno"><a href="1,55" id="srcline55"> 55</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,56" id="srcline56"> 56</a></span><span class="line">    function J = maxEntropyEnhance(I, isBad)</span></span>
<span class="srcline"><span class="lineno"><a href="1,57" id="srcline57"> 57</a></span><span class="line">        Y = rgb2gm(real(max(imresize(I, [50 50]), 0))); <span class="comment">% max - avoid complex number </span></span></span>
<span class="srcline"><span class="lineno"><a href="1,58" id="srcline58"> 58</a></span><span class="line">        </span></span>
<span class="srcline"><span class="lineno"><a href="1,59" id="srcline59"> 59</a></span><span class="line">        if exist('isBad', 'var')</span></span>
<span class="srcline"><span class="lineno"><a href="1,60" id="srcline60"> 60</a></span><span class="line">            isBad = (imresize(isBad, [50 50]));</span></span>
<span class="srcline"><span class="lineno"><a href="1,61" id="srcline61"> 61</a></span><span class="line">            Y = Y(isBad);</span></span>
<span class="srcline"><span class="lineno"><a href="1,62" id="srcline62"> 62</a></span><span class="line">        end</span></span>
<span class="srcline"><span class="lineno"><a href="1,63" id="srcline63"> 63</a></span><span class="line">        </span></span>
<span class="srcline"><span class="lineno"><a href="1,64" id="srcline64"> 64</a></span><span class="line">        if isempty(Y)</span></span>
<span class="srcline"><span class="lineno"><a href="1,65" id="srcline65"> 65</a></span><span class="line">           J = I; <span class="comment">% no enhancement k = 1</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,66" id="srcline66"> 66</a></span><span class="line">           return;</span></span>
<span class="srcline"><span class="lineno"><a href="1,67" id="srcline67"> 67</a></span><span class="line">        end</span></span>
<span class="srcline"><span class="lineno"><a href="1,68" id="srcline68"> 68</a></span><span class="line">        </span></span>
<span class="srcline"><span class="lineno"><a href="1,69" id="srcline69"> 69</a></span><span class="line">        opt_k = fminbnd(@(k) ( -entropy(applyK(Y, k)) ),1, 7);</span></span>
<span class="srcline"><span class="lineno"><a href="1,70" id="srcline70"> 70</a></span><span class="line">        J = applyK(I, opt_k, a, b) - 0.01;</span></span>
<span class="srcline"><span class="lineno"><a href="1,71" id="srcline71"> 71</a></span><span class="line">        </span></span>
<span class="srcline"><span class="lineno"><a href="1,72" id="srcline72"> 72</a></span><span class="line">    end</span></span>
<span class="srcline"><span class="lineno"><a href="1,73" id="srcline73"> 73</a></span><span class="line">end</span></span>
<span class="srcline"><span class="lineno"><a href="1,74" id="srcline74"> 74</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,75" id="srcline75"> 75</a></span><span class="line">function I = rgb2gm(I)</span></span>
<span class="srcline"><span class="lineno"><a href="1,76" id="srcline76"> 76</a></span><span class="line">if size(I,3) == 3</span></span>
<span class="srcline"><span class="lineno"><a href="1,77" id="srcline77"> 77</a></span><span class="line">    I = im2double(max(0,I)); <span class="comment">% negative double --&gt; complex double</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,78" id="srcline78"> 78</a></span><span class="line">    I = ( I(:,:,1).*I(:,:,2).*I(:,:,3) ).^(1/3);</span></span>
<span class="srcline"><span class="lineno"><a href="1,79" id="srcline79"> 79</a></span><span class="line">end</span></span>
<span class="srcline"><span class="lineno"><a href="1,80" id="srcline80"> 80</a></span><span class="line">end</span></span>
<span class="srcline"><span class="lineno"><a href="1,81" id="srcline81"> 81</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,82" id="srcline82"> 82</a></span><span class="line">function J = applyK(I, k, a, b)</span></span>
<span class="srcline"><span class="lineno"><a href="1,83" id="srcline83"> 83</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,84" id="srcline84"> 84</a></span><span class="line">if ~exist( 'a', 'var' )</span></span>
<span class="srcline"><span class="lineno"><a href="1,85" id="srcline85"> 85</a></span><span class="line">    a = -0.3293;</span></span>
<span class="srcline"><span class="lineno"><a href="1,86" id="srcline86"> 86</a></span><span class="line">end</span></span>
<span class="srcline"><span class="lineno"><a href="1,87" id="srcline87"> 87</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,88" id="srcline88"> 88</a></span><span class="line">if ~exist( 'b', 'var' )</span></span>
<span class="srcline"><span class="lineno"><a href="1,89" id="srcline89"> 89</a></span><span class="line">    b = 1.1258;</span></span>
<span class="srcline"><span class="lineno"><a href="1,90" id="srcline90"> 90</a></span><span class="line">end</span></span>
<span class="srcline"><span class="lineno"><a href="1,91" id="srcline91"> 91</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,92" id="srcline92"> 92</a></span><span class="line">f = @(x)exp((1-x.^a)*b);</span></span>
<span class="srcline"><span class="lineno"><a href="1,93" id="srcline93"> 93</a></span><span class="line">beta = f(k);</span></span>
<span class="srcline"><span class="lineno"><a href="1,94" id="srcline94"> 94</a></span><span class="line">gamma = k.^a;</span></span>
<span class="srcline"><span class="lineno"><a href="1,95" id="srcline95"> 95</a></span><span class="line">J = I.^gamma.*beta;</span></span>
<span class="srcline"><span class="lineno"><a href="1,96" id="srcline96"> 96</a></span><span class="line">end</span></span>
<span class="srcline"><span class="lineno"><a href="1,97" id="srcline97"> 97</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,98" id="srcline98"> 98</a></span><span class="line">function S = tsmooth( I, lambda, sigma, sharpness)</span></span>
<span class="srcline"><span class="lineno"><a href="1,99" id="srcline99"> 99</a></span><span class="line">if ( ~exist( 'lambda', 'var' ) )</span></span>
<span class="srcline"><span class="lineno"><a href="1,100" id="srcline100">100</a></span><span class="line">    lambda = 0.01;</span></span>
<span class="srcline"><span class="lineno"><a href="1,101" id="srcline101">101</a></span><span class="line">end</span></span>
<span class="srcline"><span class="lineno"><a href="1,102" id="srcline102">102</a></span><span class="line">if ( ~exist( 'sigma', 'var' ) )</span></span>
<span class="srcline"><span class="lineno"><a href="1,103" id="srcline103">103</a></span><span class="line">    sigma = 3.0;</span></span>
<span class="srcline"><span class="lineno"><a href="1,104" id="srcline104">104</a></span><span class="line">end</span></span>
<span class="srcline"><span class="lineno"><a href="1,105" id="srcline105">105</a></span><span class="line">if ( ~exist( 'sharpness', 'var' ) )</span></span>
<span class="srcline"><span class="lineno"><a href="1,106" id="srcline106">106</a></span><span class="line">    sharpness = 0.001;</span></span>
<span class="srcline"><span class="lineno"><a href="1,107" id="srcline107">107</a></span><span class="line">end</span></span>
<span class="srcline"><span class="lineno"><a href="1,108" id="srcline108">108</a></span><span class="line">I = im2double( I );</span></span>
<span class="srcline"><span class="lineno"><a href="1,109" id="srcline109">109</a></span><span class="line">x = I;</span></span>
<span class="srcline"><span class="lineno"><a href="1,110" id="srcline110">110</a></span><span class="line">[ wx, wy ] = computeTextureWeights( x, sigma, sharpness);</span></span>
<span class="srcline"><span class="lineno"><a href="1,111" id="srcline111">111</a></span><span class="line">S = solveLinearEquation( I, wx, wy, lambda );</span></span>
<span class="srcline"><span class="lineno"><a href="1,112" id="srcline112">112</a></span><span class="line">end</span></span>
<span class="srcline"><span class="lineno"><a href="1,113" id="srcline113">113</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,114" id="srcline114">114</a></span><span class="line">function [ W_h, W_v ] = computeTextureWeights( fin, sigma, sharpness)</span></span>
<span class="srcline"><span class="lineno"><a href="1,115" id="srcline115">115</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,116" id="srcline116">116</a></span><span class="line">dt0_v = [diff(fin,1,1);fin(1,:)-fin(end,:)];</span></span>
<span class="srcline"><span class="lineno"><a href="1,117" id="srcline117">117</a></span><span class="line">dt0_h = [diff(fin,1,2)';fin(:,1)'-fin(:,end)']';</span></span>
<span class="srcline"><span class="lineno"><a href="1,118" id="srcline118">118</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,119" id="srcline119">119</a></span><span class="line">gauker_h = filter2(ones(1,sigma),dt0_h);</span></span>
<span class="srcline"><span class="lineno"><a href="1,120" id="srcline120">120</a></span><span class="line">gauker_v = filter2(ones(sigma,1),dt0_v);</span></span>
<span class="srcline"><span class="lineno"><a href="1,121" id="srcline121">121</a></span><span class="line">W_h = 1./(abs(gauker_h).*abs(dt0_h)+sharpness);</span></span>
<span class="srcline"><span class="lineno"><a href="1,122" id="srcline122">122</a></span><span class="line">W_v = 1./(abs(gauker_v).*abs(dt0_v)+sharpness);</span></span>
<span class="srcline"><span class="lineno"><a href="1,123" id="srcline123">123</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,124" id="srcline124">124</a></span><span class="line">end</span></span>
</pre>
</div>
<pre class="code">
<span class="srcline"><span class="lineno"><a href="1,125" id="srcline125">125</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,126" id="srcline126">126</a></span><span class="line"><span class="keyword">function</span> <span class="var type1" id="S91T133U554">OUT</span> = solveLinearEquation( <span class="var type1" id="S92T133U557">IN</span>, <span class="var type1" id="S93T133U558">wx</span>, <span class="var type1" id="S94T133U559">wy</span>, <span class="var type1" id="S95T2U560">lambda</span> )</span></span>
<span class="srcline"><span class="lineno"><a href="1,127" id="srcline127">127</a></span><span class="line"><span class="mxinfo " id="T2:U6">[ <span class="var type1" id="S96T2U564">r</span>, <span class="var type1" id="S97T2U565">c</span>, <span class="var type1" id="S98T2U566">ch</span> ] = size( <span class="var type1" id="S92T133U569">IN</span> )</span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,128" id="srcline128">128</a></span><span class="line"><span class="mxinfo " id="T2:U11"><span class="var type1" id="S100T2U572">k</span> = <span class="mxinfo " id="T2:U13"><span class="var type1" id="S96T2U574">r</span> * <span class="var type1" id="S97T2U575">c</span></span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,129" id="srcline129">129</a></span><span class="line"><span class="mxinfo " id="T146:U16"><span class="var type1" id="S101T146U578">dx</span> =  <span class="mxinfo " id="T146:U18"><span class="mxinfo " id="T2:U19">-<span class="var type1" id="S95T2U581">lambda</span></span> * <span class="mxinfo " id="T146:U21"><span class="var type1" id="S93T133U583">wx</span>( : )</span></span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,130" id="srcline130">130</a></span><span class="line"><span class="mxinfo " id="T146:U23"><span class="var type1" id="S102T146U587">dy</span> =  <span class="mxinfo " id="T146:U25"><span class="mxinfo " id="T2:U26">-<span class="var type1" id="S95T2U590">lambda</span></span> * <span class="mxinfo " id="T146:U28"><span class="var type1" id="S94T133U592">wy</span>( : )</span></span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,131" id="srcline131">131</a></span><span class="line"><span class="mxinfo " id="T133:U30"><span class="var type1" id="S103T133U596">tempx</span> = <span class="mxinfo " id="T133:U32">[<span class="mxinfo " id="T86:U33"><span class="var type1" id="S93T133U600">wx</span>(<span class="mxinfo " id="T87:U35">:</span>,<span class="mxinfo " id="T2:U36">end</span>)</span>,<span class="mxinfo " id="T140:U37"><span class="var type1" id="S93T133U605">wx</span>(<span class="mxinfo " id="T87:U39">:</span>,<span class="mxinfo " id="T147:U40"><span class="mxinfo " id="T2:U41">1</span>:<span class="mxinfo " id="T2:U42">end-1</span></span>)</span>]</span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,132" id="srcline132">132</a></span><span class="line"><span class="mxinfo " id="T133:U43"><span class="var type1" id="S105T133U615">tempy</span> = <span class="mxinfo " id="T133:U45">[<span class="mxinfo " id="T130:U46"><span class="var type1" id="S94T133U619">wy</span>(<span class="mxinfo " id="T2:U48">end</span>,<span class="mxinfo " id="T122:U49">:</span>)</span>;<span class="mxinfo " id="T137:U50"><span class="var type1" id="S94T133U625">wy</span>(<span class="mxinfo " id="T148:U52"><span class="mxinfo " id="T2:U53">1</span>:<span class="mxinfo " id="T2:U54">end-1</span></span>,<span class="mxinfo " id="T122:U55">:</span>)</span>]</span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,133" id="srcline133">133</a></span><span class="line"><span class="mxinfo " id="T146:U56"><span class="var type1" id="S106T146U635">dxa</span> = <span class="mxinfo " id="T146:U58"><span class="mxinfo " id="T2:U59">-<span class="var type1" id="S95T2U638">lambda</span></span> *<span class="mxinfo " id="T146:U61"><span class="var type1" id="S103T133U640">tempx</span>(:)</span></span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,134" id="srcline134">134</a></span><span class="line"><span class="mxinfo " id="T146:U63"><span class="var type1" id="S107T146U644">dya</span> = <span class="mxinfo " id="T146:U65"><span class="mxinfo " id="T2:U66">-<span class="var type1" id="S95T2U647">lambda</span></span> *<span class="mxinfo " id="T146:U68"><span class="var type1" id="S105T133U649">tempy</span>(:)</span></span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,135" id="srcline135">135</a></span><span class="line"><span class="mxinfo " id="T133:U70"><span class="var type1" id="S103T133U653">tempx</span> = <span class="mxinfo " id="T133:U72">[<span class="mxinfo " id="T86:U73"><span class="var type1" id="S93T133U657">wx</span>(<span class="mxinfo " id="T87:U75">:</span>,<span class="mxinfo " id="T2:U76">end</span>)</span>,<span class="mxinfo " id="T140:U77">zeros(<span class="var type1" id="S96T2U663">r</span>,<span class="mxinfo " id="T2:U79"><span class="var type1" id="S97T2U665">c</span>-1</span>)</span>]</span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,136" id="srcline136">136</a></span><span class="line"><span class="mxinfo " id="T133:U81"><span class="var type1" id="S105T133U669">tempy</span> = <span class="mxinfo " id="T133:U83">[<span class="mxinfo " id="T130:U84"><span class="var type1" id="S94T133U673">wy</span>(<span class="mxinfo " id="T2:U86">end</span>,<span class="mxinfo " id="T122:U87">:</span>)</span>;<span class="mxinfo " id="T137:U88">zeros(<span class="mxinfo " id="T2:U89"><span class="var type1" id="S96T2U681">r</span>-1</span>,<span class="var type1" id="S97T2U683">c</span>)</span>]</span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,137" id="srcline137">137</a></span><span class="line"><span class="mxinfo " id="T146:U92"><span class="var type1" id="S109T146U686">dxd1</span> = <span class="mxinfo " id="T146:U94"><span class="mxinfo " id="T2:U95">-<span class="var type1" id="S95T2U689">lambda</span></span> * <span class="mxinfo " id="T146:U97"><span class="var type1" id="S103T133U691">tempx</span>(:)</span></span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,138" id="srcline138">138</a></span><span class="line"><span class="mxinfo " id="T146:U99"><span class="var type1" id="S110T146U695">dyd1</span> = <span class="mxinfo " id="T146:U101"><span class="mxinfo " id="T2:U102">-<span class="var type1" id="S95T2U698">lambda</span></span> * <span class="mxinfo " id="T146:U104"><span class="var type1" id="S105T133U700">tempy</span>(:)</span></span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,139" id="srcline139">139</a></span><span class="line"><span class="mxinfo " id="T86:U106"><span class="mxinfo " id="T86:U107"><span class="var type1" id="S93T133U705">wx</span>(<span class="mxinfo " id="T87:U109">:</span>,<span class="mxinfo " id="T36:U110">end</span>)</span> = <span class="mxinfo " id="T2:U111">0</span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,140" id="srcline140">140</a></span><span class="line"><span class="mxinfo " id="T130:U112"><span class="mxinfo " id="T130:U113"><span class="var type1" id="S94T133U713">wy</span>(<span class="mxinfo " id="T36:U115">end</span>,<span class="mxinfo " id="T122:U116">:</span>)</span> = <span class="mxinfo " id="T2:U117">0</span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,141" id="srcline141">141</a></span><span class="line"><span class="mxinfo " id="T146:U118"><span class="var type1" id="S111T146U720">dxd2</span> = <span class="mxinfo " id="T146:U120"><span class="mxinfo " id="T2:U121">-<span class="var type1" id="S95T2U723">lambda</span></span> * <span class="mxinfo " id="T146:U123"><span class="var type1" id="S93T133U725">wx</span>(:)</span></span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,142" id="srcline142">142</a></span><span class="line"><span class="mxinfo " id="T146:U125"><span class="var type1" id="S112T146U729">dyd2</span> = <span class="mxinfo " id="T146:U127"><span class="mxinfo " id="T2:U128">-<span class="var type1" id="S95T2U732">lambda</span></span> * <span class="mxinfo " id="T146:U130"><span class="var type1" id="S94T133U734">wy</span>(:)</span></span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,143" id="srcline143">143</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,144" id="srcline144">144</a></span><span class="line"><span class="var type0" id="S113T0U738">Ax</span> = <span class="message error" id="M14F659C">spdiags( <span class="mxinfo " id="T149:U132">[<span class="var type1" id="S109T146U743">dxd1</span>,<span class="var type1" id="S111T146U744">dxd2</span>]</span>, <span class="mxinfo " id="T43:U135">[<span class="mxinfo " id="T2:U136"><span class="mxinfo " id="T2:U137">-<span class="var type1" id="S100T2U749">k</span></span>+<span class="var type1" id="S96T2U750">r</span></span>,<span class="mxinfo " id="T2:U140">-<span class="var type1" id="S96T2U752">r</span></span>]</span>, <span class="var type1" id="S100T2U753">k</span>, <span class="var type1" id="S100T2U754">k</span> )</span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,145" id="srcline145">145</a></span><span class="line"><span class="var type0" id="S115T0U757">Ay</span> = <span class="message error" id="M15F659C">spdiags( <span class="mxinfo " id="T149:U144">[<span class="var type1" id="S110T146U762">dyd1</span>,<span class="var type1" id="S112T146U763">dyd2</span>]</span>, <span class="mxinfo " id="T43:U147">[<span class="mxinfo " id="T2:U148"><span class="mxinfo " id="T2:U149">-<span class="var type1" id="S96T2U768">r</span></span>+<span class="mxinfo " id="T2:U151">1</span></span>,<span class="mxinfo " id="T2:U152">-<span class="mxinfo " id="T2:U153">1</span></span>]</span>, <span class="var type1" id="S100T2U772">k</span>, <span class="var type1" id="S100T2U773">k</span> )</span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,146" id="srcline146">146</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,147" id="srcline147">147</a></span><span class="line"><span class="mxinfo " id="T146:U156"><span class="var type1" id="S116T146U776">D</span> = <span class="mxinfo " id="T146:U158"><span class="mxinfo " id="T2:U159">1</span> - ( <span class="mxinfo " id="T146:U160"><span class="mxinfo " id="T146:U161"><span class="mxinfo " id="T146:U162"><span class="var type1" id="S101T146U783">dx</span> + <span class="var type1" id="S102T146U784">dy</span></span> + <span class="var type1" id="S106T146U785">dxa</span></span> + <span class="var type1" id="S107T146U786">dya</span></span>)</span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,148" id="srcline148">148</a></span><span class="line"><span class="var type0" id="S117T0U789">A</span> = (<span class="var type0" id="S113T0U794"><span class="message error" id="M16F659C">Ax</span></span>+<span class="var type0" id="S115T0U795">Ay</span>) + (<span class="var type0" id="S113T0U799">Ax</span>+<span class="var type0" id="S115T0U800">Ay</span>)' + spdiags( <span class="var type0" id="S116T0U803">D</span>, 0, <span class="var type0" id="S100T0U805">k</span>, <span class="var type0" id="S100T0U806">k</span> );</span></span>
<span class="srcline"><span class="lineno"><a href="1,149" id="srcline149">149</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,150" id="srcline150">150</a></span><span class="line"><span class="keyword">if</span> <span class="message error" id="M17F659C">exist( <span class="mxinfo " id="T16:U167"><span class="string">'ichol'</span></span>, <span class="mxinfo " id="T19:U168"><span class="string">'builtin'</span></span> )</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,151" id="srcline151">151</a></span><span class="line">    <span class="var type0" id="S119T0U815">L</span> = ichol( <span class="var type0" id="S117T0U818"><span class="message error" id="M18F659C">A</span></span>, struct( <span class="string">'michol'</span>, <span class="string">'on'</span> ) );</span></span>
<span class="srcline"><span class="lineno"><a href="1,152" id="srcline152">152</a></span><span class="line">    <span class="mxinfo " id="T133:U169"><span class="var type1" id="S91T133U825">OUT</span> = <span class="var type1" id="S92T133U826">IN</span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,153" id="srcline153">153</a></span><span class="line">    <span class="keyword">for</span> <span class="var type1" id="S122T2U829">ii</span> = <span class="mxinfo " id="T2:U173">1:<span class="var type1" id="S98T2U832">ch</span></span></span></span>
<span class="srcline"><span class="lineno"><a href="1,154" id="srcline154">154</a></span><span class="line">        <span class="mxinfo " id="T133:U175"><span class="var type1" id="S123T133U835">tin</span> = <span class="mxinfo " id="T133:U177"><span class="var type1" id="S92T133U837">IN</span>( <span class="mxinfo " id="T87:U179">:</span>, <span class="mxinfo " id="T122:U180">:</span>, <span class="var type1" id="S122T2U840">ii</span> )</span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,155" id="srcline155">155</a></span><span class="line">        [ <span class="var type0" id="S124T0U844"><span class="message error" id="M20F659C">tout</span></span>, <span class="message error" id="M21F659C">~</span> ] = pcg( <span class="var type0" id="S117T0U848"><span class="message error" id="M19F659C">A</span></span>, <span class="var type0" id="S123T0U850">tin</span>( : ), 0.1, 50, <span class="var type0" id="S119T0U854">L</span>, <span class="var type0" id="S119T0U856">L</span>' );</span></span>
<span class="srcline"><span class="lineno"><a href="1,156" id="srcline156">156</a></span><span class="line">        <span class="var type0" id="S91T0U860">OUT</span>( :, :, <span class="var type0" id="S122T0U863">ii</span> ) = reshape( <span class="var type0" id="S124T0U866"><span class="message error" id="M22F659C">tout</span></span>, <span class="var type0" id="S96T0U867">r</span>, <span class="var type0" id="S97T0U868">c</span> );</span></span>
<span class="srcline"><span class="lineno"><a href="1,157" id="srcline157">157</a></span><span class="line">    <span class="keyword">end</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,158" id="srcline158">158</a></span><span class="line"><span class="keyword">else</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,159" id="srcline159">159</a></span><span class="line">    <span class="mxinfo " id="T133:U182"><span class="var type1" id="S91T133U872">OUT</span> = <span class="var type1" id="S92T133U873">IN</span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,160" id="srcline160">160</a></span><span class="line">    <span class="keyword">for</span> <span class="var type1" id="S122T2U876">ii</span> = <span class="mxinfo " id="T2:U186">1:<span class="var type1" id="S98T2U879">ch</span></span></span></span>
<span class="srcline"><span class="lineno"><a href="1,161" id="srcline161">161</a></span><span class="line">        <span class="mxinfo " id="T133:U188"><span class="var type1" id="S123T133U882">tin</span> = <span class="mxinfo " id="T133:U190"><span class="var type1" id="S92T133U884">IN</span>( <span class="mxinfo " id="T87:U192">:</span>, <span class="mxinfo " id="T122:U193">:</span>, <span class="var type1" id="S122T2U887">ii</span> )</span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,162" id="srcline162">162</a></span><span class="line">        <span class="var type0" id="S124T0U890">tout</span> = <span class="var type0" id="S117T0U892"><span class="message error" id="M23F659C">A</span></span>\<span class="var type0" id="S123T0U894">tin</span>( : );</span></span>
<span class="srcline"><span class="lineno"><a href="1,163" id="srcline163">163</a></span><span class="line">        <span class="var type0" id="S91T0U899">OUT</span>( :, :, <span class="var type0" id="S122T0U902">ii</span> ) = reshape( <span class="var type0" id="S124T0U905"><span class="message error" id="M24F659C">tout</span></span>, <span class="var type0" id="S96T0U906">r</span>, <span class="var type0" id="S97T0U907">c</span> );</span></span>
<span class="srcline"><span class="lineno"><a href="1,164" id="srcline164">164</a></span><span class="line">    <span class="keyword">end</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,165" id="srcline165">165</a></span><span class="line"><span class="keyword">end</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,166" id="srcline166">166</a></span><span class="line"><span class="keyword">end</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,167" id="srcline167">167</a></span><span class="line"></span></span>
</pre>
</body>
</html>
